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We study the geometrically induced cohesion of ensembles of granular "u-particles" which me- 
chanically entangle through particle interpenetration. We vary the length-to- width ratio l/w of the 
u-particles and form them into free-standing vertical columns. In laboratory experiment we monitor 
the response of the columns to sinusoidal vibration (frequency /, peak acceleration T). Column col- 
lapse occurs in a characteristic time, r, which follows the relation r = / _1 exp(A/F). A resembles 
an activation energy and is maximal at intermediate l/w. Simulation reveals that optimal strength 
results from competition between packing and entanglement. 



Many living and nonliving materials are composed of 
ensembles of non-spherical particles. In the biological 
world, ant rafts [l| and eagles nests Q are held together 
through the interpenetration of elongated and bent limbs 
and branches, respectively. In soft condensed matter sys- 
tems composed of ensembles of particles, particle shape 
influences rheological and structural properties like vis- 
cosity [1] , yield stress 0, 0| and packing density 0, @ • 

In studies of granular materials (GM), collections of 
athermal particles that interact through dissipative re- 
pulsive contact forces, most attention has been devoted 
to approximately spherical (convex) particle shapes . 
Dry non-convex granular particles also interact through 
repulsive contact forces. However non-convex particles 
possess the ability to interpenetrate, and thus depending 
on shape may experience an effective cohesion through 
interpenetration. We refer to such interpenetrating en- 
sembles as "entangled" ensembles. While the static pack- 
ing of granular non-convex particles 0, [H and the rhe- 
ology of non-convex particle suspensions Q have been 
previously studied, the role of particle entanglement in 
granular dynamics and stability is unexplored. 

In this paper we report laboratory experiment, com- 
puter simulation, and theory to elucidate the role of 
particle shape in the packing and stability of a model 
non-convex granular system of u-shaped particles ("u- 
particles"). Non-convex particles in the shape of a "u" 
are one of the simplest shapes in which concavity can be 
systematically controlled by varying arm length. Similar 
to previous studies [jl-Hlj] [l2[ of dry GM we use ver- 
tical vibration to probe relaxation and rearrangement of 
the u-particle ensembles, thereby characterizing the effec- 
tive cohesiveness of the entangled ensembles. We show 
that non-convex particle shape leads to effects which are 
qualitatively different from concave dry GM and display 
phenomena similar to those observed in studies of wet 
cohesive granular materials fllj |. 



Methods- The u-particles in our experiments consisted 
of steel staples (Duo-fast; Vernon Hills, II.) of con- 
stant width, u>=1.17 cm, and variable length, I (l/w £ 
[0.02, 1.125]), see Figs. 1 and 2. The cross section of all 
particles was rectangular with thickness of 0.5 mm and 
width 1.27 mm which corresponded to a rod-like aspect 
ratio for l/w = 0.02 particles of ps 14. Particles with 
l/w = 0.02 ± 0.02, 0.13 ± 0.02, 0.15 ± 0.03 and 0.28 ± 0.04 
were cut to length by hand and the remaining sizes were 
purchased. Columns occupied a volume V = irh(d/2) 2 , 
and the bulk volume fraction was calculated as 6 = -^tt 
where M is total particle mass and p s t — 7.85 g/cm 3 is 
the density of steel. 

Sinusoidal vibration was generated by an electromag- 
netic shaker (VTS; Aurora, OH; Fig. la). The shaker 
was controlled by Lab VIEW and a Tecron 7550 power 
amplifier. Acceleration of the vibration table was mea- 
sured by an accelerometer (PCB Piezotronics; Depew, 
New York). Columns were subjected to vibration for a 
maximum of 20 minutes or until complete column col- 
lapse occurred. 

Column collapse was monitored using a digital camera 
(Point-Grey; Richmond, BC, Canada). Image capture 
was triggered at constant phase in oscillation. Images 
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FIG. 1. Experiment and simulation of u-particle ensembles, 
a) Free-standing vertical columns are formed on a vibration 
table and the collapse dynamics are observed under vertical 
vibration, b) Computer generated ensembles of u-particles. 
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were analyzed in Matlab. 

In Monte-Carlo simulation (Fig. lb) we studied the 
packing of ensembles of u-particles. U-particles were cre- 
ated from three sphero-cylinders, each of diameter D, 
placed at right angles with backbone width W = 14D, 
measured from center to center of arms, and variable barb 
length l/w 6 [0, 1.3]. We estimate the particle volume to 
be v p = it(D/2) 2 (2L + w + §(£>/2) 2 ) The cross-sectional 
particle shape differed between simulation and experi- 
ment (rectangular in experiment, circular with diameter 
D in simulation). Excluded volume, v e , was numerically 
estimated by determining the probability p of overlap 
between two randomly oriented and positioned particles 
within a fixed spherical volume V. Overlap was deter- 
mined by computing the minimum distance between all 
line segments composing the two particles; if this distance 
was less than the sphero-cylinder diameter, D, particles 
overlapped. Excluded volume is the product of overlap 
probability and sphere volume v e = pV. 

We used a packing algorithm to generate packings of 
non-overlapping u-particles. Packing proceeded in two 
steps: first particles were placed at random position and 
orientation inside a cubic volume of cross sectional area 
(52D x 52D) . If particle placement resulted in overlap the 
placement was rejected and a new position was randomly 
selected. If after 10,000 iterations a non-overlapping par- 
ticle location was not found the algorithm proceeded to 
step two. In the second step randomly selected particles 
were displaced downwards a small random direction and 
distance ^ . If the new location of the particle resulted in 
particle overlap the particle was returned to the original 
location and a new particle chosen. The algorithm was 
halted after the center of mass height of the ensemble 
reached a steady-state. The volume fraction of the simu- 
lated packing was determined by measuring the average 
height of the pile to obtain the occupied volume and then 
dividing this by the total volume of particles. 

Column formation- In experiment, collections of 
monodisperse particles with fixed l/w were formed into 
free standing cylindrical columns with column diameter, 
d=AA cm, and height, ho = 3 cm. Columns were pre- 
pared by pouring particles into a cylindrical container 
followed by a 20 s sinusoidal vibration of the base at a 
frequency, / = 30 Hz, and peak acceleration, r = 2 (in 
units of gravitational acceleration g). We confirmed that 
steady state volume fraction was reached through our 
preparation protocol in separate experiments conducted 
over a 60 s time period. 

The final volume fraction of the column following 
packing, <pf, was a monotonically decreasing function 
of l/w (Fig. 2b). The value of 0/ observed for l/w = 
0.02 particles was within the range observed in cylin- 
drical rod experiments with comparable aspect ratio 
(length/thickness « 14) which pack to 0/ = 0.28 - 0.34 
depending upon preparation [f| @, [l3| . 

Following the packing preparation we removed the con- 
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FIG. 2. Packing of u-particles. a) <j)(t) during column prepa- 
ration, b) Final packing fraction, 4>f, as a function of par- 
ticle geometry in experiment (column diameter d = 4.4 cm 
diamonds and d = 5.6 cm squares) and simulation (white 
circles). Line is theory prediction from the random contact 
model. (Inset) Excluded volume of sphero-cylinder u-particles 
from simulation. 



fining container leaving the column free-standing (Fig. 
3a). During removal of the confining cylinder the l/w = 
0.02 particles were marginally stable with partial column 
collapse occurring approximately 50% of the time, sim- 
ilar to the results reported in [6j. Spontaneous collapse 
of the l/w > 0.02 columns was rarely observed. 

Column collapse- We applied sinusoidal vibration to 
the base of the free-standing column and observed the 
collapse process. We characterized collapse dynamics 
by monitoring the centroid height, h(t), of the column 
(Fig. 3). The collapse duration decreased with increas- 
ing r (Fig. 3b) and collapse dynamics were well described 
by a phenomenological stretched exponential fit function 
= el - ^)' 3 ]. The parameter r is the characteristic 
collapse time and ft is the stretching parameter [14j . The 
stretched exponential function is frequently applied to 
the description of relaxation dynamics of disordered sys- 
tems 



Consistent with previous studies [lS El P was in the 
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FIG. 3. Collapse dynamics, a) The normalized centroid 
height, h{t) /ho (red dot) of the column during collapse. Vi- 
bration parameters are r = 2, f = 30Hz, and images are sep- 
arated by 90 oscillation periods, b) Relaxation of h(t)/ho is 
shown for l/w = 0.379 for F = 1.23, 1.48, 1.70, 1.96, 2.20, 2.53 
respectively. Fit lines are stretched exponentials given in the 
text. 
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FIG. 4. Particle shape and optimal geometric cohesion, a) 
The logarithm of the relaxation time versus inverse accelera- 
tion with exponential fit lines r = f~ 1 e A ^ F (to = 1 s). Curves 
are offset vertically for clarity. Error bars are standard devi- 
ation of 4 or greater replicates, b) A as a function of l/w. 
Dashed line indicates estimated maximum of A (see [16]). 
Error bars represent 95% confidence interval of the best fit 
lines from (a). 



range of 0.5 - 1 and decreased slightly as T increased 
but was independent of particle geometry. For fixed l/w, 
t decreased with increasing T. The logarithm of t in- 
creased with 1 /r (Fig. 4a) and was fit by an exponential 
t = / _1 e A / r with A as the single fit parameter. 

The exponential fit is indicative of an Arrhenius-like 
process observed in the relaxation of activated systems. 
The Arrheniuse process describes the escape probability 
of a thermally or mechanically activated particle from a 
potential well of depth A. In thermal systems the escape 
time is proportional to one over the Boltzmann factor 
cxp(— -^p) where E is the activation energy required to 
overcome the potential barrier. In our system thermal 
effects are negligible, and we hypothesize that the me- 
chanical excitation plays the role of a thermal energy 
source (r analogous to kT) and A is analogous to an 
energy barrier resulting from particle entanglement. 

Column collapse occurred through the separation of 
entangled particles during vibration. We therefore ex- 
pected that the hindrance of motion due to particle 
entanglement — and thus A — would increase monotoni- 
cally with the size of the concave region and thus particle 
length. Instead we found that A was a non-monotonic 
function of l/w (Fig. 4b) with A reaching a maximum 
value at intermediate l/w = 0.394±0.045 [16]. A appears 
in an exponential and thus the relaxation time for fixed 
r displays a surprising sensitivity to variation of particle 
shape (Inset Fig. 4b). We posit that the maximum in A 
is related to the statistics of particle entanglement within 
the bulk, and we next study entanglement propensity in 
theory and simulation. 

Model and simulation- In simulation we generated u- 
particle ensembles of varying l/w (Fig. lb) and measured 
the packing fraction, 4>f. We found excellent agreement 
between experiment and simulation in (f>f (Fig. 2b) . The 
scaling of 4>f with l/w can be understood through the 
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FIG. 5. Statistics of particle entanglement in simulation, a) 
The probability distribution of entanglement number, N, as 
a function of l/w. Curves are shifted vertically for clarity; 
l/w = at top and increases in increments of 0.036 down. 
Inset: red particle containing 6 entangled neighbors, b) Mean 
values for N and p measured in simulation (circles) and the 
theoretical fit (black line), c) Density of entanglements as a 
function of l/w and the theoretical fit (black line). Vertical 
dashed line and gray bar correspond to the mean and standard 
deviation of the estimated maximum of A from experiment. 

random contact model originally developed to describe 
the packing of hard rods [8[ . 

Assuming a homogeneous packing of particles, the ran- 
dom contact model relates the excluded volume, v e , occu- 
pied volume, v p , and average number of particle contacts, 
C, to the bulk packing fraction 4>f through <j)f — C^f-. 

We numerically computed v e (Inset Fig. 2b) and found 
excellent agreement between the prediction of the ran- 
dom contact model and the values of <pf measured in ex- 
periment and simulation (Fig. 2b). Fitting the random 
contact model to experimental and simulation packing 
data we obtained C = 8.74±0.04, close to the value found 
in studies of colloid and granular rod packing 0, d, [l]| . 

We hypothesized that particle entanglement within the 
column would influence the relaxation time during ver- 
tical vibration. Thus we expected that the maximum 
in A should correspond to a maximum in the density 
of particle entanglements. In simulation we defined two 
particles as entangled when the center line of one particle 
intersected the internal plane of the neighboring particle 
(Inset Fig. 5a). We measured the number of entangle- 
ments per particle, N, for each particle in simulation. 
The probability distribution function, P(N), was sensi- 
tive to l/w (Fig. 5a) with mean value (N) increasing 
monotonically with l/w (Fig. 5b). The increase was sub- 
linear indicating that (N) grew slower than that of the 
particle's convex area (I — D)(w — 2D). 
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The scaling of (TV) with l/w can be determined by 
considering the solid volume occupied by the entangled 
particles in the focal particles convex region (the con- 
vex area with infinitesimal thickness 8). Assuming a 
homogeneous packing the solid volume in this region is 
V en t — (j>f{l ~ D){w — 2D) 5. Since each entangled par- 
ticle contributes only a portion to V en t in the shape of 
an ellipse of thickness 5, on average V en t = ol{N)w5Q- 
where a > 1 accounts for the non-planar crossings 
(Inset Fig. 5a). Solving the above relations yields 

(N) = • With single fit parameter, 

a = 2.648 ±0.108, we find excellent agreement between 
the predicted number of entanglements per particle and 
those measured in simulation (Fig. 5b). 

The spatial density of particle entanglements is p ent — 
(N)p where p = — is the number density of parti- 
cles (Fig. 5b). Substitution for (N) yields p en t = 
4£_ / v p (i 2D) ^ an( j again the simulation and the- 

ory are in good agreement (Fig. 5c) using the previosly 
determined fit parameters C and a. Furthermore the ex- 
perimental maximum max[A] at l/w — 0.394 ± 0.045 is 
close to the value obtained in simulation and theory of 
l/w = 0.340 ± 0.015 suggesting that the large relaxation 
times for the intermediate u-particle columns is due to 
the large density of mechanical entanglements. 

Conclusion- Similar to rod-like particles 0, d, [l3j |. 
columns formed from u-particles are stabilized through 
the inhibition of particle rotation and translation among 
the entangled particles. The addition of the transverse 
ends which form concave u-particles leads to mechanical 
entanglement and increases column stability. However 
the increase in entanglement with increasing length is 
offset by the decrease in particle packing density. These 
two trends conspire to generate a maximum in the density 
of mechanical entanglements in collections of non-convex 
particles of intermediate I / w-thus columns of these par- 
ticles most strongly resist separation. 

The random contact model utilized to explain the op- 
timum geometry for entanglement of u-shaped particles 
assumes only uncorrelated particle contacts within the 
bulk. Thus we expect the results to apply to rigid non- 
convex particulate systems of all scales. For instance a 
recent study of suspension rheology found that convex 
particles of differing shape collapsed to a viscosity-stress 
master curve while concave particle did not collapse to 
this curve; this difference was attributed to particle en- 
tanglement effects Q. At the micro-scale, polymers with 
rigid pendants oriented perpendicular to the polymer 
chain increase internal molecular free volume and hinder 
polymer motion which significantly affects rheology sim- 



proof rafts and shelters-which have been described as 
akin to living chain mail (l9| -through the interlocking 
and entanglement of limbs and mandibles [1]. 
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ilar to geometric entanglement UJ. At the macroscale 
strain-stiffe ning of model polymers is associated with en- 
tanglement [l8j. Even organisms can benefit from geo- 
metric entanglement. For example the fire ant Solenopsis 
invicta and the army ant Eciton burchelli create water- 
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